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Abstract 

We study a discrete two-dimensional nonlinear system that allows 
for discrete breather solutions. We perform a spectral analysis of the 
lattice dynamics at thermal equilibrium and use a cooling technique to 
measure the amount of breathers at thermal equilibrium. Our results 
confirm the existence of an energy threshold for discrete breathers. 
The cooling method provides with a novel computational technique of 
measuring and analyzing discrete breather distribution properties in 
thermal equilibrium. 

1 Introduction 

Nonlinear discrete systems support discrete breathers (DBs). These time- 
periodic and spatially localized solutions are the result of the interplay be- 
tween nonlinearity and discreteness jTj. Many studies of DBs have been 
successfully launched, on such topics as rigorous existence proofs, dynam- 
ical and structural stability and computational methods of obtaining DBs 
in classical models as well as their quantum aspects. In addition DBs have 
been detected and studied experimentally in such different systems as in- 
teracting Josephson junction systems 01, coupled nonlinear optical waveg- 
uides 3 , lattice vibrations in crystals [I], antiferromagnetic structures jS], 
micromechanical cantilever arrays 6 , Bose-Einstein condensates loaded on 
optical lattices [7], layered high-T c superconductors [S]. Their existence is 
also predicted to exist in the dynamics of dusty plasma crystals jH] • 

Among several intriguing unresolved questions concerning DBs, a cen- 
tral issue is the contribution of DBs to the dynamics of systems at thermal 
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equilibrium. Indeed breather-like excitations have been observed in a vari- 
ety of different models at finite temperatures ^U]. For some special mod- 
els with additional conservation laws semianalytical statements about the 
contribution of DBs to thermal equilibrium have been derived The 
question is then whether we can identify the contribution of DBs to vari- 
ous equilibrium and relaxation properties (like e.g. the charge trapping in 
DNA |12j). Since DBs are dynamical excitations, their contribution will be 
observable mainly in time-dependent (or frequency-dependent) correlation 
functions, while static correlation functions, e.g. specific heat, only probe 
the available energy landscape, and are not suitable for detecting dynamical 
correlations. A good way to proceed is to use a specific property of DBs and 
trace its contribution to correlation functions. Another question is whether 
we can design computational methods to separate the breather excitations 
at a given time from the rest of the excitations in the lattice. This would 
allow us to perform systematic studies of distribution properties of DBs in 
a given lattice at a given temperature. 

Concerning the first point from above, we know that in general nonlinear 
systems in two and three dimensions support breather solutions that have a 
positive lower energy threshold This is a very important property that 
can be of help in the detection of DBs in experiments. In fact for specific 
one-dimensional lattices this property holds as well. However the search for 
traces of these thresholds in correlation functions for one-dimensional lattices 
at thermal equilibrium turned out to be very complicated ^3]. There are 
two reasons for that. First, DBs in one-dimensional lattices act as strong 
scatterers of plane waves |15) . Consequently radiation can be efficiently 
trapped between DBs, and contribute to a strong interaction between DBs. 
It is also hard to find a way for letting the radiation out of the system (see 
|16j). Secondly frequency-dependent correlation functions probe gaps in 
frequency space. Although the existence of energy thresholds for breathers 
leads also to frequency thresholds, the values of these frequency thresholds 
may become too small to be easily detected [Ti] . 

In this work we examine a two-dimensional lattice. First fingerprints 
of the DB energy threshold in two-dimensional lattices have been reported 
in J2j and also in ^Sj. A comprehensive study of DBs in thermal equilib- 
rium and their influence on various relaxation mechanisms was provided by 
Ivanchenko et al 1!) . We attempt to go sufficiently beyond these studies. 
Here we present not only the frequency dependence of correlation functions, 
but also their temperature (or energy) dependence. We then proceed to ap- 
ply a cooling technique at the boundaries of our system to efficiently get rid 
of extended excitations in the lattice. The remaining localized breather-like 
excitations can then be easily analyzed. We observe the existence of energy 
thresholds and provide with novel distribution functions for DB energies. 
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2 Gap determination in the Hamiltonian system 



We study a two-dimensional quadratic lattice with one degree of freedom 
per site. The equation of motion for the particle at site in a lattice of 
size N x N with free ends is given by 

ilij = k(u i+ ij + Ui-ij - 2uij) + k(uij + i + Uij-x - 2uij) - - ufj, (1) 

where Uij is the displacement of the particle at site The masses of 
the particles are equal to unity and the interaction between them is har- 
monic with strength k while the on-site potential is of hard-(/> 4 type. The 
corresponding Hamiltonian from which Q is derived, is given by 

1 1 1 k 

H = Y. hi d > hi >J = ^Ij + 2 U i,j + 1 u tj + 4 ~ Ul ^ 2 ■ ( 2 ) 

i,j NN 

Here hij is the discrete on-site energy density which will be used later, and 
NN stands for nearest neighbors and implies (/ — i) 2 + (m — j) 2 = 1. For 
small amplitudes the linearized equation yields plane waves ui m (t) ~ 
e i(u q t-q x l-q y m) the phonon dispersion relation 

4=l + 4^(sin 2 (|) + sin 2 (|)) (3) 

where q = (q x , q y ) is the wave vector of the reciprocal lattice. 

We obtain the energy threshold for breathers using the Newton algo- 
rithm and the anticontinuous limit pp. We compute breather solutions for 
different frequencies for coupling k = 0.05. The phonon band of small am- 
plitude plane waves is located between the frequencies Wg- = ( .o) = = 1 and 
u^ = t^^\ = uj n = 1.183. We use a lattice of 19 x 19 sites. In Fig. ^we plot 
the total energy of the breather solution as a function of its frequency Qf,. 
Additionally we check the stability of the obtained solutions using the stan- 
dard Floquet analysis pQ. We observe in Fig. [J that the curve of the total 
energy of the breather as a function of frequency consists of two branches, 
separated by a minimum at frequency Slb.th- I n the left branch (dashed line) 
DBs are unstable according to the Floquet analysis (see Fig|5J^) while in 
the right branch (solid line) DBs are linearly stable solutions (FigEP). The 
breather with frequency Qb,th = 1-207 is separating the frequency region 
lu-x < f2& < £lb,th where breathers are unstable and the region VL^th < ^f> 
where breathers are linearly stable according to our computations. This spe- 
cific breather solution has a total energy E^th = 0.268. 60% of the energy of 
this DB is located on the central site, i.e. h centra i = 0.167. In other words, 
the minimum energy breather is still a rather discrete object, involving es- 
sentially only a few lattice sites in its dynamics (Fig|2J$) (see also jl3|.|20|). 
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Figure 1: Energy of the breather solution for coupling k = 0.05 as a function 
of its frequency (circles). Lines connecting circles are guides to the eye. The 
dashed vertical lines indicate the boundaries of the phonon band. 
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Figure 2: Displacements U[ m (t = 0) of breather solutions at initial time 
with zero velocities ui m (t = 0) = 0. Below each profile the eigenvalues of 
the Floquet matrix of linearized perturbations are shown in the complex 
plane (squares) together with the reference unit circle. A: = 1.188; B: 
n b = 1.207; C: Q b = 1.319. 

In terms of energy the threshold value E b ^h provides with a rigorous 
lower bound - no breathers exist with energies less than E^th- Consequently 
we may expect this feature to be detectable in the temperature dependence 
of equilibrium correlation functions. The corresponding order of magnitude 
of the average energy per site is expected to be in the region of values 
0.05—0.2. While the upper value is simply close to E bj th, the lower one can be 
obtained by observing that essentially one central cite and four neighboring 
sites are important for the dynamics of stable breathers, providing with 
a lower estimated value of E^h/b. At the same time the frequency gap 
a;,,- < Qb < Q b ,th provides with a less rigorous bound. We can only speculate 
that breathers with frequencies belonging to this gap region are less probable 
to be excited, because they are linearly unstable. Nevertheless they do exist, 
and in thermal equilibrium the system may be excited close to these solutions 
for short times. 

3 Thermal equilibrium 

We thermalize the lattice using as a first method Langevin equations of 
motion . We add to the right hand side of equations Q damping terms 
—jiii.j and a Gaussian white noise force £(t). The friction is chosen to be 
7 = 0.01. The Gaussian white noise £ is characterized by the standard 
correlation function (£(£)£(£')) = 2 r ykBT5(t — t') where ks = 1 is the dimen- 
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Figure 3: Dependence of the average energy per site {h[ m ) on the tempera- 
ture in the parameter range of interest. 

sionless Boltzmann constant and T is the temperature. We simulate the 
Langevin equations until the averaged kinetic energy per particle is close to 
the desired value T/2. We then switch off the friction and Gaussian white 
noise terms, and continue integrating the Hamiltonian equations ([T]). Time 
is set to t = at this moment of switching from Langevin to Hamiltonian 
evolution. We then measure the time and ensemble averaged kinetic energy 
and reobtain the corresponding value of T, together with computing the 
time and ensemble averaged energy per site. The temperature values are 
very close to the corresponding averaged total energies per site (/t/, m ) as 
shown in Fig|21 

We compute the time dependent displacement-displacement correlation 
function, 

11 1 N N rtfin 

= -^wh ^EE / Mt + t'^iAt'Vdt'}, (4) 

£>0 ^ tfi n -Un^p^ J 

where the brackets denote an ensemble average. In this study the ensemble 
averaging is obtained by averaging the results for 10 different realizations 
or runs. The parameters are ti n = 8000 and tfi n = 8500. The Fourier 
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Figure 4: Power spectra P(uj) as a function of frequency for various values 
of temperature. The coupling between the sites is k = 0.05. The values of 
temperature are indicated in the labels in each subgraph. The two vertical 
dashed lines mark the band edge values uiq and uj n of the phonon band while 
the vertical solid line marks Qbth = 1-207. 



transform of S(t) is: 

jfin y 

S(u m ) = At e~^ mAtJ S(Atj), (5) 

3=0 

where uj m = f 20 ^, m = 0,1, , h^^?- — 1. Finally the power spectrum of 

''fin * 

the correlation function S(t) is given by 

P{u) = {Re{S{u))f + {Im{S{u))f, (6) 

where Re(S(u>)) and Im(S(u))) denote the real and imaginary part respec- 
tively. The parameter So is chosen such that S(t = 0) = 1 implying that 
the frequency integral over S(u>) is constant. 

In Fig El we plot the power spectrum P(u) for various temperatures. 
For low temperature values the power spectrum shows that the modes with 
frequencies inside the phonon band are excited. As the temperature in- 
creases the power spectrum is shifted to larger frequencies due to the hard 
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anharmonicity. Additionally we observe that for temperatures T ~ 0.04 
and higher the spectrum exhibits a different behavior in the frequency gap 
region, viz. shows a decrement near the gap. Let us examine the power 
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Figure 5: Power spectra P(T) as a function of temperature for 4 different 
frequencies. From top to the bottom: frequency inside the phonon band, 
u) = 1.156; frequency inside the frequency gap where breathers are unstable 
to = 1.194; frequency lies outside the two previous regions to = 1.257; largest 
frequency value to = 1.508. 

spectra in detail. In Fig. [5] we plot the value of the power spectrum P(T) 
as a function of temperature for four specific frequencies. The top panel of 
Fig. |S] is plotted for a frequency inside the phonon band, the second panel 
is plotted for a frequency that belongs to the breather frequency gap, while 
the lower two panels are plotted for frequencies that are located outside the 
phonon band and the frequency gap region. In the top panel of the Fig. 
Elwe observe a maximum as the temperature increases due to the shift of 
the peak of the phonon spectrum with temperature in FigEJ This peak is 
located at u = 1.156 for low temperatures. With increasing temperature the 
peak passes through the specific frequency that we chose in the top panel of 
FigEl The characteristic temperature value of this peak will depend on the 
chosen value of the probe frequency inside the phonon band. 

In the second panel in FigEl where the chosen frequency belongs to the 
breather frequency gap, a crossover is observed at T ~ 0.1. This crossover is 
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even more evident in the third panel. In this case the frequency lies outside 
the gap region but, at the same time, is very close to it. A crossover is 
observed again at T ~ 0.1, a value that has the same order of magnitude as 
the energy threshold value E^fa that we found in Section 2. Finally in the 
bottom panel in FigEJa trace of this threshold is still observed, although the 
probe frequency is located far from the phonon band and breather frequency 
gap values, and the power spectrum P{oS) has only small tails there. 

Thus we do observe a characteristic crossover feature in the power spec- 
tra, but not as initially expected upon varying the frequency at a fixed 
temperature. Instead we observe the crossover for a fixed probe frequency 
by varying the temperature. The absence of a clearly observable frequency 
gap may be due to the fact that the frequency gap itself is rather narrow, 
that unstable breathers may still be excited and persist for some time with 
some probability, and that anharmonic extended waves contribute as well. 
An interesting question is why the observed crossover upon varying the tem- 
perature is showing a maximum or similar cusp for frequencies close to the 
frequency gap, and an (although bearly visible) opposite behavior at larger 
probe frequencies (as seen in the bottom panel in FigJEJ. A possible answer 
is that for low temperatures only occasionally breathers are excited, and the 
shift of frequency contributions to higher values with increasing tempera- 
ture is as well caused by simple hard anharmonicity effects. However once 
the temperature reaches the threshold value for breathers, breathers with 
larger energies are easily excited as well. This will cause a depletion of the 
power spectra at lower frequencies, at the expense of the increase at higher 
frequencies, just as observed in Fig03 Still the above argumentation shows 
that it is hard to separate the contribution of anharmonic phonons from 
that of breathers. We clearly are in need of a technique which does this 
separation. The next chapter will provide with a solution to this problem. 

4 How to measure breathers in thermal equilib- 
rium 

As already discussed in the introduction, DBs act as strong point-like scat- 
terers of plane waves |15| . In one-dimensional lattices this circumstance 
makes it hard to let the delocalized radiative part of excitations out of the 
system |19| . However in two-dimensional (and even more efficiently in three- 
dimensional) systems point-like scatterers will not hinder plane waves from 
moving around such an obstacle. Consequently we may attach dissipative 
boundaries to our lattice, and expect the radiation to disappear from the 
system during a reasonable short time, leaving the immobile localized excita- 
tions behind. Such cooling techniques have been used for studying breather 
relaxations |16j , JH] > [El ■ Here we are only interested in letting the radiation 
out and being left with DBs. 
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We study an ensemble of 10 different realizations that correspond to the 
same initial energy per site Eq. We thermalize the system in another way 
here. We start the system with a given initial energy Eq, i.e. we integrate 
the system for 1000 time units using initially random displacements for 
the particles. After that dissipative boundaries are switched on for a given 
transient of time. This is done here by adding a friction term to the boundary 
sites of a 20 x 20 lattice with friction constant 0.1 for some cooling time T CO oi- 
While much more sophisticated ways of reflectionless dissipative boundaries 
can be implemented (U, this simple method suffices for the results presented 
below. 

The effect of the dissipative boundary is shown in Figs l6l7l where the 
energy per site remaining in the lattice is plotted as a function of time, both 
in normal and in logarithmic units. Note that the energy values are normal- 
ized to their corresponding number at time t = at which the dissipative 
boundary is switched on. 
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Figure 6: The time dependence of the normalized lattice energy per site 
in the presence of a dissipative boundary for different values of the initial 
energy per site as indicated in the figures. 

We clearly observe that after some transient behavior the lattice energy 
density changes rather slowly, indicating the expected outcome - delocalized 
excitations left the system, leaving the localized ones behind. Both FigEl 
and especially Fig|3show that the characteristic waiting time increases with 
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increasing initial energy density from roughly T coo \ = 2000 up to T coo [ = 
10000. This can be observed from the shift of the inflection points in FigEJ 
In the following we will use a cooling time T coo i = 10 4 , but we will also 
compare with shorter cooling times. 

In Fig. El we plot the total energy that remains in the lattice (excluding 
the sites that belong to the boundaries) per site and per number of realiza- 
tions E as a function of initial energy Eq. Note the log scale for the y-axis. 
For clarity the inset is for the same data but with linear axis scaling. We 
observe a crossover for Eq ~ 0.05. We fit the curve of E(Eq) first using 

E b,th 

E = AE b , th e E o (7) 
with parameters A and E^^. This holds assuming that only breathers with 

E b,th 

E(,th win be excited. Then the probability to form a DB is e E o while its 

E b,th 

contribution to an energy distribution will be E^ ^e E o . Depending on 
the energy range we use for fitting we obtain Ej,^ = 0.17 (Eq < 0.06) or 
Eb,th = 0.233 (Eq < 0.1). We were not able to fit data at larger energies. 
The obtained values 0.17-0.233 are in good agreement with the expected 
value of the DB energy threshold 0.27. The dashed line in FigEl is the 
corresponding fit with parameters A = 0.175 and E^ ^ = 0.198. 
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Figure 8: Graph of the average energy per site and per number of realizations 
E as a function of the initial energy Eq after cooling. The error bars are 
the standard deviations of the mean for the statistical ensemble. Note the 
logarithmic scale of the y-axis. Solid line -guide to the eye. Dashed line - 
fit using ((7J). Dashed-dotted line - fit using (jSJ) (for parameters see text). 
Inset: Same but with linear y-axis scaling to observe the crossover around 
E = 0.05. 

In a refined fitting we allow also for breathers with larger energies. Then 
the energy per site contribution to E will be given by 



where p{E^) is the density of DB states. Assuming p(Eb) = 1 (its actual con- 

E b,th 

stant value can be always absorbed in A), Q yields A(EoEb :t h + E^)e E » 
and the subsequent fitting procedure results similar results for E^^ as ob- 
tained with ®. The dashed-dotted line in FiglHl is the corresponding fit 
with parameters A = 5.53 and E^^ = 0.163. While the low energy region 
of the numerical data is well reproduced, the high energy data are underes- 
timated with (J2J) and overestimated with (jHJ). This implies that the density 
p(Eb) is not constant, but decaying with increasing energy E^. 

But we can obtain even more relevant data. After the cooling process we 
measure all values h^ m of the energy density. In Figs. |§] and EH we present 




(8) 
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their distribution W{E) additionally averaged over all the realizations after 
cooling. The initial energy per site that is used is indicated in the labels 
in each subfigure. Note that the data are obtained by coarsegraining the 
energy axis with a grid size of 0.05. The data obtained within the first box 
< E < 0.05 are omitted from the plots (except for the two left upper 
panels at the two lowest energies Eq). 

The DB energy threshold is visible in all the subgraphs except the first 
two subgraphs where the initial energy Eq is very low and DBs cannot be 
formed. We also clearly see that already for energies Eq = 0.04 a whole 
distribution of breather energies is obtained, with maximum breather en- 
ergies being twice larger compared to the minimum DB energy. Note that 
the observed energy gaps are very close to the value 0.17. This corresponds 
exactly to the expected energy of the central site for the minimum energy 
breather being equal to 0.167 (see Fig|2j3). The nearest neighbours of such 
a breather carry approximately ten percent of the full DB energy each, i.e. 
0.027. Their contributions are thus within the mentioned first box and not 
present in the plotted data. 
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Figure 9: Graph of the energy distribution W{E). See text for details. 
The arrows have length 0.17 and show the energy region below which the 
distribution is expected to vanish due to the DB energy threshold. 
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In order to test the sensitivity of these data on the cooling time T coo \ we 
present in Figs 1 1 1 1 f21 similar results for T coo i = 5000 which is twice smaller. 
While the overall statistics needs to be improved in both cases, we find semi- 
quantitative agreement. This indicates that during the rather long cooling 
times the statistical properties of the localized excitations remaining in the 
lattice are not significantly changing. The observed increase of W(E) for 
small E and large temperatures may be due to the fact that at these large 
temperatures the cooling time T coo \ was too short, leaving delocalized ex- 
citations inside the system. We will provide with further evidence for the 
correctness of this conclusion. 

The above analysis suggested that to some accuracy the excitations in 
the lattice at thermal equilibrium can be considered as a sum of localized 
and delocalized excitations. The power spectra in Fig0]represent thus a sum 
of the power spectra of both types of excitations. We use now the cooling 
process, which leaves us with the localized excitations only. We wait for 
T CO oi = 8000 and after that compute the power spectrum of the remaining 
localized excitations in the system. The results are shown in Figs ll3|lH 
The comparison of Figs ll3lHl with Fig0] shows, that for energies Eq < 0.03 
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i) the dynamics is essentially governed by delocalized excitations; ii) after 
the cooling period only delocalized states with almost zero group velocities 
remain in the system (these correspond to the band edge frequencies and to 
the peak inside the band). For larger energies Eq the power spectra after 
cooling show the existence of localized excitations with frequencies outside 
the band uj'i. Moerover, we observe very clearly that the frequency gap 
< f) 6 < &bth is depleted here, indicating that unstable DBs decay during 
the cooling time and radiate their energy into the dissipative boundary, 
even if they were present to some extent in thermal equilibrium. Thus the 
presented way of separating localized from delocalized excitations in thermal 
equilibrium allows also to quantitatively determine the contribution of DBs 
to correlation functions. Finally we observe that at the largest temperatures 
considered, some nonzero statistical weight is observed in the frequency 
region of the phonon band and even below it. This confirms our previous 
expectation, that at these large temperatures the cooling times are probably 
not long enough to let the extended excitations out of the system. A possible 
reason may be that at these large temperatures many DBs are coexisting in 
the system. Extended waves will spend more time to diffuse around these 
scattering centers. DBs could even form temporary percolation networks 
which would hinder the escape of waves even more efficiently. 



5 Conclusions 

Our results show that discrete breathers (DBs) leave clear and detectable 
fingerprints in the thermal equilibrium properties of nonlinear lattices. Us- 
ing the case of a two-dimensional lattice, we demonstrated the persistence 
of an energy threshold for the existence of DBs in time-dependent corre- 
lation functions by observing weak crossover features. Moreover, we used 
the technique of boundary cooling to separate the localized excitation part 
at thermal equilibrium from the delocalized ones. This allows us to study 
statistical properties of DB excitations in thermal equilibrium, e.g. their 
contribution to the abovementioned correlation functions. This simple step 
allowed us to unambiguously confirm the presence of the energy threshold 
for DBs in thermal equilibrium. We furthermore confirm that the frequency 
gap (which corresponds to the excitation of unstable DBs) is depleted after 
the separation procedure. Thus we conclude that unstable DBs are not con- 
tributing to dynamical correlations at thermal equilibrium. By comparing 
the spectra of the system at thermal equilibrium with the spectra of the 
DB part only, we can provide with reliable statistical weights of both the 
delocalized and localized excitation parts of the system at thermal equilib- 
rium at one and the same frequency. It seems to be very difficult to provide 
with a similar cooling and separation procedure in any realistic experimental 
setup, despite the fact that DBs have been observed in a variety of different 
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systems. That makes computational studies of DB properties at thermal 
equilibrium a unique way of gaining further understanding of the relaxation 
and excitation properties of complex lattices. 

With this simple technique, which certainly needs more refinement, we 
are now able to reliably compute various statistical contributions of DBs, 
including also correlation effects between DBs and delocalized excitations. 
The door is thus open for starting serious analytical work on DBs in thermal 
equilibrium, since the presented new numerical tools will allow for a much 
more refined testing of various theories as compared to simulations of ther- 
mal equilibrium only. 
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